#ifndef _MGFACTORY_H_
#define _MGFACTORY_H_

#include "Multigrid/Multigrid.h"
#include <string>
#include <fstream>
#include <jsoncpp/json/json.h>

template <int Dim, class RST, class ITPL> class MGFactory;

template <int Dim, class RST, class ITPL>
class MGFactory{
public:

  using iVec = Vec<int,2>;

  using rVec = Vec<Real,2>;

  template <class T>
  using Vector = std::vector<T>;

  MGFactory(const RectDomain<Dim>& aDomain, const char* aBCTypes);

  MGFactory(const Vec<Real,Dim> &alo, const Vec<Real,Dim> &ahi, const
            int an, const char* aBCTypes);

  MGFactory(const std::string& jsonfile);
  
  void setParam(const Real aweight, const int anumIter1, const int
                anumIter2, const int anumIter3, const Real areltol,
                const int amaxIter);

  const RectDomain<Dim>& getDomain() const;

  const char* getBCTypes() const;

  const Real getWeight() const;

  const Vec<int,3> getNumIter() const;

  const Real getReltol() const;

  const int getMaxIter() const;
  
  template <class TFunc>
  Tensor<Real,Dim> SolveVCycle(const TFunc& func);

  template <class TFunc>
  Tensor<Real,Dim> SolveFMVCycle(const TFunc& func);

  template <class TFunc>
  Real computeError(const Tensor<Real,Dim>& res, const TFunc& func,
                    const int p = 0);

  void plot(const Tensor<Real,Dim>& res, const std::string &file);
protected:
  Multigrid<Dim> MG; 
};

template <int Dim, class RST, class ITPL>
MGFactory<Dim,RST,ITPL>::MGFactory(const RectDomain<Dim>& aDomain, const char* aBCTypes){
  Vector<RectDomain<Dim> > vD{aDomain};
  int an = (aDomain.size())[0]-1;
  for(int n = an ; n >= 4 ; n/=2){
    RectDomain<Dim> Domain = vD.back().coarsen();
    vD.push_back(Domain);
  }
  int size = vD.size();
  Vector<Restriction<Dim>*> vR;
  Vector<Interpolation<Dim>*> vI;
  for (int i = 0 ; i < size - 1 ; i++){
    Restriction<Dim>* Rst = new RST(vD[i],vD[i+1]);
    vR.push_back(Rst);
    Interpolation<Dim>* Itpl = new ITPL(vD[i],vD[i+1]);
    vI.push_back(Itpl);
  }
  MG = Multigrid<Dim>(vD,vR,vI);
  MG.setBCType(aBCTypes);
}

template <int Dim, class RST, class ITPL>
MGFactory<Dim,RST,ITPL>::MGFactory(const Vec<Real,Dim> &alo, const Vec<Real,Dim> &ahi, const
                                   int an, const char* aBCTypes){
  Vec<int,Dim> lo(0);
  Vec<int,Dim> hi(an);
  Box<Dim> bx(lo,hi);
  Vec<Real,Dim> dx = (ahi-alo)/an;
  RectDomain<Dim> aDomain(bx,dx,NodeCentered,1);
  Vector<RectDomain<Dim> > vD{aDomain};
  for(int n = an ; n >= 4 ; n/=2){
    RectDomain<Dim> Domain = vD.back().coarsen();
    vD.push_back(Domain);
  }
  int size = vD.size();
  Vector<Restriction<Dim>*> vR;
  Vector<Interpolation<Dim>*> vI;
  for (int i = 0 ; i < size - 1 ; i++){
    Restriction<Dim>* Rst = new RST(vD[i],vD[i+1]);
    vR.push_back(Rst);
    Interpolation<Dim>* Itpl = new ITPL(vD[i],vD[i+1]);
    vI.push_back(Itpl);
  }
  MG = Multigrid<Dim>(vD,vR,vI);
  MG.setBCType(aBCTypes);
}

template <int Dim, class RST, class ITPL>
MGFactory<Dim,RST,ITPL>::MGFactory(const std::string& jsonfile){
  Json::Value root;
  Json::Reader reader;
  std::ifstream os;
  os.open(jsonfile);
  reader.parse(os, root);
  os.close();
  Json::Value Domain = root["Domain"];
  Vec<Real,Dim> rlo;
  Vec<Real,Dim> rhi;
  switch(Dim){
  case 1:
    rlo = Vec<Real,Dim>{Domain["lo"][0].asDouble()};
    rhi = Vec<Real,Dim>{Domain["hi"][0].asDouble()};
    break;
  case 2:
    rlo = Vec<Real,Dim>{Domain["lo"][0].asDouble(),Domain["lo"][1].asDouble()};
    rhi = Vec<Real,Dim>{Domain["hi"][0].asDouble(),Domain["hi"][1].asDouble()};
    break;
  }
  int an = Domain["n"].asInt();
  const char* bctypes = root["bctypes"].asString().c_str();
  char BCTypes[4];
  switch(Dim){
  case 1:
    BCTypes[0] = bctypes[0];BCTypes[1] = bctypes[1];
    break;
  case 2:
    BCTypes[0] = bctypes[0];BCTypes[1] = bctypes[1];
    BCTypes[2] = bctypes[2];BCTypes[3] = bctypes[3];
    break;
  }
  //std::string restriction = root["restriction"].asString();
  //std::string interpolation = root["interpolation"].asString();
  const Real weight = root["weight"].asDouble();
  const int numiter1 = root["numiter1"].asInt();
  const int numiter2 = root["numiter2"].asInt();
  const int numiter3 = root["numiter3"].asInt();
  const int maxiter = root["maxiter"].asInt();
  const Real reltol = root["reltol"].asDouble();
  Vec<int,Dim> lo(0);
  Vec<int,Dim> hi(an);
  Box<Dim> bx(lo,hi);
  Vec<Real,Dim> dx = (rhi-rlo)/an;
  RectDomain<Dim> aDomain(bx,dx,NodeCentered,1);
  Vector<RectDomain<Dim> > vD{aDomain};
  for(int n = an ; n >= 4 ; n/=2){
    RectDomain<Dim> Domain = vD.back().coarsen();
    vD.push_back(Domain);
  }
  int size = vD.size();
  Vector<Restriction<Dim>*> vR;
  Vector<Interpolation<Dim>*> vI;
  for (int i = 0 ; i < size - 1 ; i++){
    Restriction<Dim>* Rst = new RST(vD[i],vD[i+1]);
    vR.push_back(Rst);
    Interpolation<Dim>* Itpl = new ITPL(vD[i],vD[i+1]);
    vI.push_back(Itpl);
  }
  MG = Multigrid<Dim>(vD,vR,vI);
  MG.setBCType(BCTypes);
  MG.setParam(weight,numiter1,numiter2,numiter3,reltol,maxiter);
}

template <int Dim, class RST, class ITPL>
void MGFactory<Dim,RST,ITPL>::setParam(const Real aweight, const int anumIter1, const int
                                       anumIter2, const int anumIter3, const Real areltol,
                                       const int amaxIter){
  MG.setParam(aweight,anumIter1,anumIter2,anumIter3,areltol,amaxIter);
}

template <int Dim, class RST, class ITPL>
const RectDomain<Dim>& MGFactory<Dim,RST,ITPL>::getDomain() const{
  return (MG.getvDomain())[0];
}

template <int Dim, class RST, class ITPL>
const char* MGFactory<Dim,RST,ITPL>::getBCTypes() const{
  return MG.getBCTypes();
}

template <int Dim, class RST, class ITPL>
const Real MGFactory<Dim,RST,ITPL>::getWeight() const{
  return MG.getWeight();
}

template <int Dim, class RST, class ITPL>
const Vec<int,3> MGFactory<Dim,RST,ITPL>::getNumIter() const{
  return MG.getNumIter();
}

template <int Dim, class RST, class ITPL>
const Real MGFactory<Dim,RST,ITPL>::getReltol() const{
  return MG.getReltol();
}

template <int Dim, class RST, class ITPL>
const int MGFactory<Dim,RST,ITPL>::getMaxIter() const{
  return MG.getMaxIter();
}

template <int Dim, class RST, class ITPL>
template <class TFunc>
Tensor<Real,Dim> MGFactory<Dim,RST,ITPL>::SolveVCycle(const TFunc& func){
  const RectDomain<Dim>& Domain = (MG.getvDomain())[0];
  FuncFiller<Dim> FuncF(Domain);
  GhostFiller<Dim> GhostF(Domain);
  Tensor<Real,Dim> rhs(Domain.getGhostedBox());
  FuncF.fillFunc(rhs,func);
  const char* BCTypes = MG.getBCTypes();
  GhostF.fillAllSides(rhs,BCTypes,func);
  Tensor<Real,Dim> phi(Domain.getGhostedBox());
  MG.SolveVCycle(phi,rhs,func);
  return phi;
}

template <int Dim, class RST, class ITPL>
template <class TFunc>
Tensor<Real,Dim> MGFactory<Dim,RST,ITPL>::SolveFMVCycle(const TFunc& func){
  const RectDomain<Dim>& Domain = (MG.getvDomain())[0];
  FuncFiller<Dim> FuncF(Domain);
  GhostFiller<Dim> GhostF(Domain);
  Tensor<Real,Dim> rhs(Domain.getGhostedBox());
  FuncF.fillFunc(rhs,func);
  const char* BCTypes = MG.getBCTypes();
  GhostF.fillAllSides(rhs,BCTypes,func);
  Tensor<Real,Dim> phi(Domain.getGhostedBox());
  MG.SolveFMVCycle(phi,rhs,func);
  return phi;
}

template <int Dim, class RST, class ITPL>
template <class TFunc>
Real MGFactory<Dim,RST,ITPL>::computeError(const Tensor<Real,Dim>& res,
                                           const TFunc& func, const int
                                           p){
  const RectDomain<Dim>& Domain = (MG.getvDomain())[0];
  Tensor<Real,Dim> exactres(Domain);
  const Vec<int,Dim>& lo = Domain.lo();
  const Vec<Real,Dim>& dx = Domain.spacing();
  switch(Dim){
  case 1:
    loop_box_1(Domain,i){
      Vec<int,Dim> Node{i};
      Vec<Real,Dim> rNode = (Node-lo)*dx;
      exactres(i) = func.F(rNode);
    }
    break;
  case 2:
    loop_box_2(Domain,i,j){
      Vec<int,Dim> Node{i,j};
      Vec<Real,Dim> rNode = (Node-lo)*dx;
      exactres(i,j) = func.F(rNode);
    }
    break;
  }
  Tensor<Real,Dim> Rres = res.slice(Domain);
  PossionOP<Dim> POP(Domain);
  return POP.computeNorm(Rres-exactres,p);
}


template <int Dim, class RST, class ITPL>
void MGFactory<Dim,RST,ITPL>::plot(const Tensor<Real,Dim>& res, const std::string &file){
  const RectDomain<Dim>& Domain = (MG.getvDomain())[0];
  Tensor<Real,Dim> Rres = res.slice(Domain);
  const Vec<int,Dim> &lo = Domain.lo();
  const Vec<Real,Dim> &dx = Domain.spacing();
  std::ofstream os;
  os.open(file);
  os << "res=[\n" << std::endl;
  switch(Dim){
  case 1:
    loop_box_1(Domain,i){
      Vec<int,Dim> Node{i};
      Vec<Real,Dim> rNode = (Node - lo)*dx;
      os << rNode[0] << "," << Rres(i) << ";\n" \
         << std::endl;
    }
    os << "];\n" << std::endl;
    os << "x=res(:,1);y=res(:,2);\n" << std::endl;
    os << "plot(x,y);\n" << std::endl;
    break;
  case 2:
    loop_box_2(Domain,i,j){
      Vec<int,Dim> Node{i,j};
      Vec<Real,Dim> rNode = (Node - lo)*dx;
      os << rNode[0] << "," << rNode[1] << "," << Rres(i,j) << ";\n"    \
         << std::endl;
    }
    os << "];\n" << std::endl;
    os << "x=res(:,1);y=res(:,2);z=res(:,3);\n" << std::endl;
    os << "T = delaunay(x,y);trimesh(T,x,y,z);\n" << std::endl;
    break;
  }
  os.close();
}


#endif // _MGFACTORY_H_
